###################################
###################################
## Analysis for Replication, JEPP submission ###
###################################
###################################

rm(list = ls())

# Define a vector with all the required package names
packages <- c(
  "rstatix", "ggpubr", "sjPlot", "interactions", "stargazer",
  "effects", "survey", "collapse", "stats", "scales",
  "dplyr", "ggthemes", "patchwork", "cjoint"
)

# Install any packages that are not already installed
install_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if(length(install_packages)) install.packages(install_packages)

# Load all the required packages
lapply(packages, library, character.only = TRUE)


library(rstatix)
library(ggpubr)
library(sjPlot)
library(interactions)
library(stargazer)
library(effects) 
library(survey)
library(collapse) 
library(stats) 
library(scales)
library(dplyr)
library(ggthemes)
library(patchwork)
library(cjoint)

#### Set working directory here 

load("all_merged.RData")
head(all) 
nrow(all)


table(all$country)

### to compare support for quota vs parity, restrict to treated

T1 <- subset(all, (Treatment== "Quota" | Treatment== "Parity"))
nrow(T1)
is.factor(T1$Treatment)
table(T1$country)


##############################
##############################
## Figure 2 ##################
##############################
##############################

### use survey to find the weighted means 

T1$fSupport<-as.factor(T1$Support)
T1$fTreatment<-as.factor(T1$Treatment)

## subset to UK to get weighted mean levels of support 
T1UK <- subset(T1, (country== "UK"))
nrow(T1UK)

### overall support 

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-T1UK[,vars]
head(T2)
wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)
svymean(~fTreatment, wdat, na.rm = TRUE)

###

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1UK[,vars], fTreatment=="Parity")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)

#         mean     SE
#fSupport1 0.15500 0.0185
#fSupport2 0.24224 0.0205
#fSupport3 0.21443 0.0201
#fSupport4 0.24718 0.0198
#fSupport5 0.14115 0.0170

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1UK[,vars], fTreatment=="Quota")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)

#           mean     SE
#fSupport1 0.13509 0.0156
#fSupport2 0.22579 0.0193
#fSupport3 0.24317 0.0212
#fSupport4 0.27136 0.0205
#fSupport5 0.12459 0.0145

cur_df <- T1UK %>%
group_by(Treatment, Support) %>% # NB: the order of the grouping
summarise(count=n()) %>% mutate(perc=count/sum(count))

cur_df<-na.omit(cur_df)
cur_df

cur_df$wperc<-NA
cur_df$wperc[cur_df$Treatment=="Parity" & cur_df$Support==1]<-0.15500
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==2]<-0.24224
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==3]<-0.21443
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==4]<-0.24718
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==5]<-0.14115
cur_df$wperc[cur_df$Treatment=="Quota" & cur_df$Support==1]<-0.13509 
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==2]<-0.22579
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==3]<-0.24317
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==4]<-0.27136
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==5]<-0.12459
cur_df

cur_df$wperc<-cur_df$wperc*100

p<- ggplot(data=cur_df, aes(x=factor(Support), y=wperc, fill=Treatment)) +
geom_bar(width=0.7, position=position_dodge(width=0.75), stat="identity") +
scale_y_continuous(labels = scales::percent_format(accuracy = 5L)) +
labs(y = "Percentage", x="To what extent do you support or oppose 
gender quota [parity] laws in your country?", title="UK") +
theme_bw() + 
 ylim(0, 40)+ 
theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank(), plot.title = element_text(hjust = 0.5)) +
theme(legend.title=element_blank())

p2<- p + scale_x_discrete(breaks=c("1", "2", "3", "4", "5"),
        labels=c("Strongly 
oppose", "Tend to 
oppose", "Neither support 
nor oppose", "Tend to 
support", "Strongly 
support"))
p2

####### Do it in greyscale

p <- ggplot(data = cur_df, aes(x = factor(Support), y = wperc, fill = Treatment)) +
  geom_bar(width = 0.7, position = position_dodge(width = 0.75), stat = "identity", color = "black") +  # Removed fill parameter
  scale_fill_manual(values = c("darkgrey", "lightgrey")) +  # Set fill colors manually
  scale_y_continuous(labels = scales::percent_format(accuracy = 5L)) +
  labs(y = "Percentage", x = "To what extent do you support or oppose gender quota [parity] laws in your country?", title = "UK") +
  theme_bw() +
  ylim(0, 40) +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(hjust = 0.5),
    legend.title = element_blank()
  )

p2<- p + scale_x_discrete(breaks=c("1", "2", "3", "4", "5"),
        labels=c("Strongly 
oppose", "Tend to 
oppose", "Neither support 
nor oppose", "Tend to 
support", "Strongly 
support"))
p2


########### Now the same figure for France 


## subset to FR to get weighted mean levels of support 
T1FR <- subset(T1, (country== "FR"))
nrow(T1FR)

### overall support 

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-T1FR[,vars]
head(T2)
wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)
svymean(~fTreatment, wdat, na.rm = TRUE)


### now for parity 

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1FR[,vars], fTreatment=="Parity")
head(T2)
wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)

#              mean     SE
#fSupport1 0.037614 0.0116
#fSupport2 0.102189 0.0217
#fSupport3 0.271401 0.0275
#fSupport4 0.357146 0.0322
#fSupport5 0.231650 0.0257


vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1FR[,vars], fTreatment=="Quota")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)

#            mean     SE
#fSupport1 0.075917 0.0191
#fSupport2 0.134288 0.0246
#fSupport3 0.234174 0.0260
#fSupport4 0.369977 0.0316
#fSupport5 0.185643 0.0253


cur_df <- T1FR %>%
group_by(Treatment, Support) %>% # NB: the order of the grouping
summarise(count=n()) %>% mutate(perc=count/sum(count))

cur_df<-na.omit(cur_df)
cur_df

cur_df$wperc<-NA
cur_df$wperc[cur_df$Treatment=="Parity" & cur_df$Support==1]<-0.037614
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==2]<-0.102189
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==3]<-0.271401
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==4]<-0.357146
cur_df$wperc[cur_df$Treatment=="Parity"& cur_df$Support==5]<-0.231650
cur_df$wperc[cur_df$Treatment=="Quota" & cur_df$Support==1]<-0.075917
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==2]<-0.134288
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==3]<-0.234174
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==4]<-0.369977
cur_df$wperc[cur_df$Treatment=="Quota"& cur_df$Support==5]<-0.185643
cur_df

cur_df$wperc<-cur_df$wperc*100

p<- ggplot(data=cur_df, aes(x=factor(Support), y=wperc, fill=Treatment)) +
geom_bar(width=0.7, position=position_dodge(width=0.75), stat="identity") +
scale_y_continuous(labels = scales::percent_format(accuracy = 5L)) +
labs(y = "Percentage", x="To what extent do you support or oppose 
gender quota [parity] laws in your country?", title="France") +
theme_bw()  +
ylim(0, 40)+ 
theme(panel.grid.minor.x=element_blank(), panel.grid.major.x=element_blank(), plot.title = element_text(hjust = 0.5)) +
theme(legend.title=element_blank())

p2<- p + scale_x_discrete(breaks=c("1", "2", "3", "4", "5"),
        labels=c("Strongly 
oppose", "Tend to 
oppose", "Neither support 
nor oppose", "Tend to 
support", "Strongly 
support"))
p2

##### do it in greyscale

p <- ggplot(data = cur_df, aes(x = factor(Support), y = wperc, fill = Treatment)) +
  geom_bar(width = 0.7, position = position_dodge(width = 0.75), stat = "identity", color = "black") +  # Removed fill parameter
  scale_fill_manual(values = c("darkgrey", "lightgrey")) +  # Set fill colors manually
  scale_y_continuous(labels = scales::percent_format(accuracy = 5L)) +
  labs(y = "Percentage", x = "To what extent do you support or oppose gender quota [parity] laws in your country?", title = "FR") +
  theme_bw() +
  ylim(0, 40) +
  theme(
    panel.grid.minor.x = element_blank(),
    panel.grid.major.x = element_blank(),
    plot.title = element_text(hjust = 0.5),
    legend.title = element_blank()
  )

p2<- p + scale_x_discrete(breaks=c("1", "2", "3", "4", "5"),
        labels=c("Strongly 
oppose", "Tend to 
oppose", "Neither support 
nor oppose", "Tend to 
support", "Strongly 
support"))
p2


##############################
##############################
## H1, Figure 1###############
##############################
##############################

## subset to UK to get weighted mean levels of support 
T1UK <- subset(T1, (country== "UK"))
nrow(T1UK)

xm <- weighted.mean(T1UK$Support[T1UK$Treatment=="Parity"], T1UK$Weight[T1UK$Treatment=="Parity"], na.rm=TRUE)
xm ###2.98
xmb <- weighted.mean(T1UK$Support[T1UK$Treatment=="Quota"], T1UK$Weight[T1UK$Treatment=="Quota"], na.rm=TRUE)
xmb ## 3.02

res1 <- lm(Support ~ Treatment, data = subset(T1, country=="UK"), weights = Weight)
summary(res1)
plot_model(res1, type = "pred")
##Predicted values (marginal effects) for specific model terms

eff <- allEffects(res1)
eff_df <- as.data.frame(eff[["Treatment"]])

##

p <- ggplot(eff_df, aes(x = Treatment, y = fit, ymin = lower, ymax = upper, linetype = Treatment)) +
  geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4, show.legend = FALSE) +  # Adding points
  labs(x = "", y = "Support (1-5)") +
  coord_flip() +
  theme_classic() +
  ylim(1, 5) +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16), guide = "none") +  # Remove legend for shape
  ggtitle("UK") +
  theme_bw()
p

p<- p +  theme(legend.position = "none", 
                 axis.text = element_text(size = 18), 
                 axis.text.x = element_text(size = 18),  # Adjust x-axis label size
                 plot.title = element_text(size = 18))
p

## subset to FR to get weighted mean levels of support 
T1FR <- subset(T1, (country== "FR"))
nrow(T1FR)

xm2 <- weighted.mean(T1FR$Support[T1FR$Treatment=="Parity"], T1FR$Weight[T1FR$Treatment=="Parity"], na.rm=TRUE)
xm2 ### 3.64
xm2b <- weighted.mean(T1FR$Support[T1FR$Treatment=="Quota"], T1FR$Weight[T1FR$Treatment=="Quota"], na.rm=TRUE)
xm2b ## 3.46

res1b <- lm(Support ~ Treatment, data = subset(T1, country=="FR"), weights = Weight)
summary(res1b)
plot_model(res1b, type = "pred")

eff <- allEffects(res1b)
eff_df <- as.data.frame(eff[["Treatment"]])

p1 <- ggplot(eff_df, aes(x = Treatment, y = fit, ymin = lower, ymax = upper, linetype = Treatment)) +
  geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4, show.legend = FALSE) +  # Adding points
  labs(x = "", y = "Support (1-5)") +
  coord_flip() +
  theme_classic() +
  ylim(1, 5) +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16), guide = "none") +  # Remove legend for shape
  ggtitle("France") +
  theme_bw()

##p1<- p1 + theme(legend.position = "none", axis.text = element_text(size = 18))
p1 <- p1 + theme(legend.position = "none", 
                 axis.text = element_text(size = 18), 
                 axis.text.x = element_text(size = 18),  # Adjust x-axis label size
                 plot.title = element_text(size = 18))

library(cowplot)
combined_plot <- plot_grid(p, p1,  ncol = 1)
combined_plot

##############################
##############################
## H2 ########################
##############################
##############################

res <- t.test(Support ~ country, data=T1)
res ## support for PA measures is significantly higher in FR compared to UK 

##############################
##############################
## H3, Figure 3 ##############
##############################
##############################

#### First UK 

### mean levels of support for men and women 

xm <- weighted.mean(T1UK$Support[T1UK$resp_gender =="Male"], T1UK$Weight[T1UK$resp_gender =="Male"], na.rm=TRUE)
xm ##   2.844172
xm <- weighted.mean(T1UK$Support[T1UK$Treatment=="Parity" & T1UK$resp_gender =="Male"], T1UK$Weight[T1UK$Treatment=="Parity"& T1UK$resp_gender =="Male"], na.rm=TRUE)
xm ##  2.779216
xm <- weighted.mean(T1UK$Support[T1UK$Treatment=="Quota"& T1UK$resp_gender =="Male"], T1UK$Weight[T1UK$Treatment=="Quota"& T1UK$resp_gender =="Male"], na.rm=TRUE)
xm ## 2.90358
xm <- weighted.mean(T1UK$Support[T1UK$resp_gender =="Female"], T1UK$Weight[T1UK$resp_gender =="Female"], na.rm=TRUE)
xm ##   3.154259
xm <- weighted.mean(T1UK$Support[T1UK$Treatment=="Parity" & T1UK$resp_gender =="Female"], T1UK$Weight[T1UK$Treatment=="Parity"& T1UK$resp_gender =="Female"], na.rm=TRUE)
xm ## 3.159968
xm <- weighted.mean(T1UK$Support[T1UK$Treatment=="Quota"& T1UK$resp_gender =="Female"], T1UK$Weight[T1UK$Treatment=="Quota"& T1UK$resp_gender =="Female"], na.rm=TRUE)
xm ##  3.148465

## overall levels - what percent of men and women support 
T1UKMEN <- subset(T1, (country== "UK" & resp_gender=="Male"))
nrow(T1UKMEN)
T1UKWOMEN <- subset(T1, (country== "UK" & resp_gender=="Female"))
nrow(T1UKWOMEN)

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1UKMEN[,vars], fTreatment=="Parity")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)
##22.3 + 10.7

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1UKMEN[,vars], fTreatment=="Quota")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1UKWOMEN[,vars], fTreatment=="Parity")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)
##22.3 + 10.7

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1UKWOMEN[,vars], fTreatment=="Quota")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)


res5 <- lm(Support ~ resp_gender, data = T1UK, weights = Weight)
summary(res5)

res5 <- lm(Support ~ Treatment + resp_gender + resp_gender:Treatment, data = T1UK, weights = Weight)
summary(res5)

eff <- allEffects(res5)
eff_df <- as.data.frame(eff[["Treatment:resp_gender"]])

p <- ggplot(eff_df, aes(x = Treatment, y = fit, color=resp_gender, ymin = lower, ymax = upper, linetype = Treatment)) +
     geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4) +  # Adding points
    labs(x = "", y = "") +
    coord_flip() +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16), guide = "none") +  # Remove legend for shape
  scale_color_manual(values = c("black", "gray")) +
    theme_classic() +
    ylim(1,5)+
  theme_bw()
    
p<- p + guides(color=guide_legend(title="Respondent Gender")) + theme(legend.position = c(0.85, 0.87),
legend.background = element_rect(fill = "white", color = "black"), , axis.text=element_text(size=18)) + 
ggtitle("UK")
p

### Then France 

xm <- weighted.mean(T1FR$Support[T1FR$resp_gender =="Male"], T1FR$Weight[T1FR$resp_gender =="Male"], na.rm=TRUE)
xm ##  3.412486
xm <- weighted.mean(T1FR$Support[T1FR$Treatment=="Parity" & T1FR$resp_gender =="Male"], T1FR$Weight[T1FR$Treatment=="Parity"& T1FR$resp_gender =="Male"], na.rm=TRUE)
xm ## 3.509119
xm <- weighted.mean(T1FR$Support[T1FR$Treatment=="Quota"& T1FR$resp_gender =="Male"], T1FR$Weight[T1FR$Treatment=="Quota"& T1FR$resp_gender =="Male"], na.rm=TRUE)
xm ## 3.316076
xm <- weighted.mean(T1FR$Support[T1FR$resp_gender =="Female"], T1FR$Weight[T1FR$resp_gender =="Female"], na.rm=TRUE)
xm ##  3.697407
xm <- weighted.mean(T1FR$Support[T1FR$Treatment=="Parity" & T1FR$resp_gender =="Female"], T1FR$Weight[T1FR$Treatment=="Parity"& T1FR$resp_gender =="Female"], na.rm=TRUE)
xm ## 3.776185
xm <- weighted.mean(T1FR$Support[T1FR$Treatment=="Quota"& T1FR$resp_gender =="Female"], T1FR$Weight[T1FR$Treatment=="Quota"& T1FR$resp_gender =="Female"], na.rm=TRUE)
xm ##  3.609612


## overall levels - what percent of men and women support 
T1FRMEN <- subset(T1, (country== "FR" & resp_gender=="Male"))
nrow(T1FRMEN)
T1FRWOMEN <- subset(T1, (country== "FR" & resp_gender=="Female"))
nrow(T1FRWOMEN)

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1FRMEN[,vars], fTreatment=="Parity")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)
##22.3 + 10.7

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1FRMEN[,vars], fTreatment=="Quota")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1FRWOMEN[,vars], fTreatment=="Parity")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)
##22.3 + 10.7

vars<-c("ID", "Weight", "fSupport", "fTreatment")
T2<-subset(T1FRWOMEN[,vars], fTreatment=="Quota")
head(T2)

wdat <- svydesign(id      = ~ID,
                          weights = ~Weight,
                          nest    = TRUE,
                          data    = T2)
svymean(~fSupport, wdat, na.rm = TRUE)


res5 <- lm(Support ~ resp_gender, data = T1FR, weights = Weight)
summary(res5)

res5 <- lm(Support ~ Treatment + resp_gender + resp_gender:Treatment, data = T1FR, weights = Weight)
summary(res5)

eff <- allEffects(res5)
eff_df <- as.data.frame(eff[["Treatment:resp_gender"]])

p2 <- ggplot(eff_df, aes(x = Treatment, y = fit, color=resp_gender, ymin = lower, ymax = upper, linetype = Treatment)) +
     geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4) +  # Adding points
    labs(x = "", y = "") +
    coord_flip() +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16), guide = "none") +  # Remove legend for shape
  scale_color_manual(values = c("black", "gray")) +
    theme_classic() +
    ylim(1,5)+
  theme_bw()
    
p2<- p2 + guides(color=guide_legend(title="Respondent Gender")) + theme(legend.position = c(0.85, 0.87),
legend.background = element_rect(fill = "white", color = "black"), , axis.text=element_text(size=18)) + 
ggtitle("France")
p2

combined_plot2 <-  plot_grid(p, p2,  ncol = 1)

combined_plot2

##############################
##############################
## H4, Figure 4 ##############
##############################
##############################
library(ggplot2)
library(sjPlot)

# UK plot
res13 <- lm(Support ~ Treatment + resp_gender + Masculine + Treatment*resp_gender*Masculine, data = T1UK, weights = Weight)
p1 <- plot_model(res13, type = "pred", terms = c("Masculine", "Treatment", "resp_gender"), colors = "bw", title = "UK", legend.title = "Treatment") +
  theme_classic()

# Modify line thickness and confidence interval fill color for UK plot
p1 <- p1 +
  scale_color_manual(values = c("Quota" = "black", "Parity" = "black")) +  # Set line color to black for both treatments
  scale_fill_manual(values = c("Quota" = "#CCCCCC", "Parity" = "#333333")) +  # Set confidence interval fill colors
  geom_line(size = 2)  # Increase line thickness

# France plot
res14 <- lm(Support ~ Treatment + resp_gender + Masculine + Treatment*resp_gender*Masculine, data = T1FR, weights = Weight)
p2 <- plot_model(res14, type = "pred", terms = c("Masculine", "Treatment", "resp_gender"), colors = "bw", title = "France", legend.title = "Treatment") +
  theme_classic()

p2 <- p2 +
  scale_color_manual(values = c("Quota" = "black", "Parity" = "black")) +  # Set line color to black for both treatments
  scale_fill_manual(values = c("Quota" = "#CCCCCC", "Parity" = "#333333")) +  # Set confidence interval fill colors
  geom_line(size = 2)  # Increase line thickness

# Combine plots
combined_plot2 <- plot_grid(p1, p2, ncol = 1)
combined_plot2

##############################
##############################
## H5, Figure 5 ##############
##############################
##############################

## open the conjoint data 
load("yougov_uk_conjoint.RData")
head(long) 
longuk<-long
nrow(longuk) #3468

res1 <- lm(Qualified_score ~ Treatment, data = subset(longuk, feature1=="Female"), weights = Weight)
summary(res1)

eff <- allEffects(res1)
eff_df <- as.data.frame(eff[["Treatment"]])

p <- ggplot(eff_df, aes(x = Treatment, y = fit, ymin = lower, ymax = upper, linetype = Treatment)) +
  geom_linerange(size = 1.5, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4, show.legend = FALSE) +  # Adding points
  labs(x = "", y = "Perceived `Qualified to be an MP' (1-11) by Treatment, Women Candidates") +
  coord_flip() +
  theme_classic() +
  ylim(3, 9) +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  ggtitle("UK") +
  theme_bw()

##p1<- p1 + theme(legend.position = "none", axis.text = element_text(size = 18))
p <- p + theme(legend.position = "none", 
                 axis.text = element_text(size = 18), 
                 axis.text.x = element_text(size = 18),  # Adjust x-axis label size
                 plot.title = element_text(size = 18))
p

## now the french conjoint data 

load("FRCONJdataset.RData")
head(long)
longfr<-long

res2 <- lm(Qualified_score ~ Treatment, data = subset(longfr, feature1=="Female"), weights = Weight)
summary(res2)

eff <- allEffects(res2)
eff_df <- as.data.frame(eff[["Treatment"]])

p2 <- ggplot(eff_df, aes(x = Treatment, y = fit, ymin = lower, ymax = upper, linetype = Treatment)) +
  geom_linerange(size = 1.5, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4, show.legend = FALSE) +  # Adding points
  labs(x = "", y = "Perceived `Qualified to be an MP' (1-11) by Treatment, Women Candidates") +
  coord_flip() +
  theme_classic() +
  ylim(3, 9) +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  ggtitle("FR") +
  theme_bw()

##p1<- p1 + theme(legend.position = "none", axis.text = element_text(size = 18))
p2 <- p2 + theme(legend.position = "none", 
                 axis.text = element_text(size = 18), 
                 axis.text.x = element_text(size = 18),  # Adjust x-axis label size
                 plot.title = element_text(size = 18))
p2

#### Figure 5
combined_plot3 <-  plot_grid(p, p2, ncol=1)
combined_plot3

###################
### Now a figure for men candidates for Appendix (Fig S12) 
###################

res1a <- lm(Qualified_score ~ Treatment, data = subset(longuk, feature1=="Male"), weights = Weight)
summary(res1a)

eff <- allEffects(res1a)
eff_df <- as.data.frame(eff[["Treatment"]])

pa <- ggplot(eff_df, aes(x = Treatment, y = fit, ymin = lower, ymax = upper, linetype = Treatment)) +
  geom_linerange(size = 1.5, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4, show.legend = FALSE) +  # Adding points
  labs(x = "", y = "Perceived `Qualified to be an MP' (1-11) by Treatment, Men Candidates") +
  coord_flip() +
  theme_classic() +
  ylim(3, 9) +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  ggtitle("UK") +
  theme_bw()

pa <- pa + theme(legend.position = "none", 
                 axis.text = element_text(size = 18), 
                 axis.text.x = element_text(size = 18),  # Adjust x-axis label size
                 plot.title = element_text(size = 18))
pa


res2a <- lm(Qualified_score ~ Treatment, data = subset(longfr, feature1=="Male"), weights = Weight)
summary(res2a)

eff <- allEffects(res2a)
eff_df <- as.data.frame(eff[["Treatment"]])

p2a <- ggplot(eff_df, aes(x = Treatment, y = fit, ymin = lower, ymax = upper, linetype = Treatment)) +
  geom_linerange(size = 1.5, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4, show.legend = FALSE) +  # Adding points
  labs(x = "", y = "Perceived `Qualified to be an MP' (1-11) by Treatment, Men Candidates") +
  coord_flip() +
  theme_classic() +
  ylim(3, 9) +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  ggtitle("FR") +
  theme_bw()

p2a <- p2a + theme(legend.position = "none", 
                 axis.text = element_text(size = 18), 
                 axis.text.x = element_text(size = 18),  # Adjust x-axis label size
                 plot.title = element_text(size = 18))

combined_plot3a <- plot_grid(pa, p2a, ncol=1)
combined_plot3a

##############################
##############################
## H6, Figure 6 ##############
##############################
##############################

res5 <- lm(Qualified_score ~ Treatment + resp_gender + resp_gender:Treatment, data = subset(longuk, feature1=="Female"), weights= Weight)
summary(res5)

eff <- allEffects(res5)
eff_df <- as.data.frame(eff[["Treatment:resp_gender"]])


p5 <- ggplot(eff_df, aes(x = Treatment, y = fit, color=resp_gender, ymin = lower, ymax = upper, linetype = Treatment)) +
     geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4) +  # Adding points
  labs(x = "", y = "Perceived `Qualified to be an MP' (1-11) by Treatment, Women Candidates") +
    coord_flip() +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  scale_color_manual(values = c("black", "gray")) +
    theme_classic() +
    ylim(3,9)+
  theme_bw()
    
p5<- p5 + guides(color=guide_legend(title="Respondent Gender")) + theme(legend.position = c(0.93, 0.87),
legend.background = element_rect(fill = "white", color = "black"), , axis.text=element_text(size=18)) + 
ggtitle("UK")
p5


res5a <- lm(Qualified_score ~ Treatment + resp_gender + resp_gender:Treatment, data = subset(longfr, feature1=="Female"), weights= Weight)
summary(res5a)

eff <- allEffects(res5a)
eff_df <- as.data.frame(eff[["Treatment:resp_gender"]])

p5a <- ggplot(eff_df, aes(x = Treatment, y = fit, color=resp_gender, ymin = lower, ymax = upper, linetype = Treatment)) +
     geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4) +  # Adding points
  labs(x = "", y = "Perceived `Qualified to be an MP' (1-11) by Treatment, Women Candidates") +
    coord_flip() +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  scale_color_manual(values = c("black", "gray")) +
    theme_classic() +
    ylim(3,9)+
  theme_bw()
    
p5a<- p5a + guides(color=guide_legend(title="Respondent Gender")) + theme(legend.position = c(0.93, 0.87),
legend.background = element_rect(fill = "white", color = "black"), , axis.text=element_text(size=18)) + 
ggtitle("FR")
p5a



combined_plot4 <- plot_grid(p5, p5a, ncol=1)
combined_plot4

##########################
##########################
### List experiment, Figure 7 ######
##########################
##########################

subset_fr <- all[all$country == "FR", ]
subset_uk <- all[all$country == "UK", ]


# UK

subset_uk$list_treated2 <- ifelse(subset_uk$list_treated == 0, "Control", "Treated")
subset_uk$list_treated2 <- factor(subset_uk$list_treated2)

# FR

subset_fr$list_treated2 <- ifelse(subset_fr$list_treated == 0, "Control", "Treated")
subset_fr$list_treated2 <- factor(subset_fr$list_treated2)

# UK
subset_uk$list_treated2<-as.factor(subset_uk$list_treated2)
subset_uk$Treatment2<-as.factor(subset_uk$Treatment2)
subset_uk$Treatment<-as.factor(subset_uk$Treatment)

res9 <- lm(lcount ~ list_treated2 + Treatment + list_treated2:Treatment, weight=Weight, data=subset_uk)
summary(res9)

eff <- allEffects(res9)
eff_df <- as.data.frame(eff[["list_treated2:Treatment"]])

p6 <- ggplot(eff_df, aes(x = Treatment, y = fit, color=list_treated2, ymin = lower, ymax = upper, linetype = Treatment)) +
     geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4) +  # Adding points
  labs(x = "", y = "") +
    coord_flip() +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  scale_color_manual(values = c("black", "gray")) +
    theme_classic() +
    ylim(1,5)+
  theme_bw()
    
p6<- p6 + guides(color=guide_legend(title="Treated (List Experiment)")) + theme(legend.position = c(0.1, 0.87),
legend.background = element_rect(fill = "white", color = "black"), , axis.text=element_text(size=18)) + 
ggtitle("UK")
p6


# FR

subset_fr$list_treated2<-as.factor(subset_fr$list_treated2)
subset_fr$Treatment2<-as.factor(subset_fr$Treatment2)
subset_fr$Treatment<-as.factor(subset_fr$Treatment)

res10 <- lm(lcount ~ list_treated2 + Treatment + list_treated2:Treatment, weight=Weight, data=subset_fr)
summary(res10)

eff <- allEffects(res10)
eff_df <- as.data.frame(eff[["list_treated2:Treatment"]])

p <- ggplot(eff_df, aes(x = Treatment, y = fit, color=list_treated2)) +
    geom_point() + 
    geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.1, size=1) +
    labs(x = "", y = "") +
    coord_flip() +
  scale_color_manual(values = c("black", "gray")) +
    theme_classic() +
    ylim(1,5)+
  theme_bw()
    
p7<- p + guides(color=guide_legend(title="Treated (List Experiment)")) + theme(legend.position = c(0.1, 0.87),
legend.background = element_rect(fill = "white", color = "black"), , axis.text=element_text(size=12)) + 
ggtitle("FR")

p7 <- ggplot(eff_df, aes(x = Treatment, y = fit, color=list_treated2, ymin = lower, ymax = upper, linetype = Treatment)) +
     geom_linerange(size = 2, show.legend = FALSE) +
  geom_point(aes(shape = Treatment), size = 4) +  # Adding points
  labs(x = "", y = "") +
    coord_flip() +
  scale_linetype_manual(values = c("Quota" = "dotted", "Parity" = "solid", "Control" = "dashed"), guide = "none") +  # Remove legend for linetype
  scale_shape_manual(values = c("Quota" = 16, "Parity" = 16, "Control" = 16), guide = "none") +  # Remove legend for shape
  scale_color_manual(values = c("black", "gray")) +
    theme_classic() +
    ylim(1,5)+
  theme_bw()
    
p7<- p7 + guides(color=guide_legend(title="Treated (List Experiment)")) + theme(legend.position = c(0.1, 0.87),
legend.background = element_rect(fill = "white", color = "black"), , axis.text=element_text(size=18)) + 
ggtitle("FR")
p7

# no differences by framing for either country
## Figure 7 

combined_plot7 <- plot_grid(p6, p7, ncol=1)
combined_plot7

